Interaction matrix element fluctuations in ballistic quantum dots: dynamical effects 
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We study matrix element fluctuations of the two-body screened Coulomb interaction and of the 
one-body surface charge potential in ballistic quantum dots, comparing behavior in actual chaotic 
billiards with analytic results previously obtained in a normalized random wave model. We find 
that the matrix element variances in actual chaotic billiards typically exceed by a factor of 3 or 
4 the predictions of the random wave model, for dot sizes commonly used in experiments. We 
discuss dynamical effects that are responsible for this enhancement. These dynamical effects have 
an even more striking effect on the covariance, which changes sign when compared with random 
wave predictions. In billiards that do not display hard chaos, an even larger enhancement of matrix 
element fluctuations is possible. These enhanced fluctuations have implications for peak spacing 
statistics and spectral scrambling for quantum dots in the Coulomb blockade regime. 
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I. INTRODUCTION 

The statistical fluctuations of single-particle energies 
and wave functions of dots whose single-particle dynam- 
ics are chaotic may be well approximated by random ma- 
trix theory (RMT) 1 . The mesoscopic fluctuations of the 
conductance through open dots that are strongly coupled 
to leads are then successfully described by RMT 2 . In the 
opposite limit of an almost-isolated dot, the charge is 
quantized and electron-electron interactions modify the 
mesoscopic fluctuations of the conductance. 

The randomness of the single-particle wave functions 
induces randomness into the interaction matrix elements 
when the latter are expressed in the basis of the former. 
These matrix elements can be decomposed into an aver- 
age and a fluctuating part. The average part of the inter- 
action, when combined with the one-body kinetic energy 
and a confining potential, leads to the so-called universal 
Hamiltonian 3 ' 4 . This universal Hamiltonian includes a 
charging energy term and an exchange interaction term 
that is proportional to the square of the total spin of the 
dot (an additional Cooper-channel term is repulsive in a 
quantum dot and can be ignored). The fluctuating part 
of the interaction is suppressed by the Thouless conduc- 
tance <7r, and in the limit gx — > oo, the dot is completely 
described by the universal Hamiltonian. 

The charging energy term leads to charge quantiza- 
tion in a weakly coupled dot, and the conductance peak 
height distributions in such a dot were derived in Ref. HI 
using the RMT statistics of the single-particle wave func- 
tions. Qualitative features of these peak height distribu- 
tions as well as the parametric peak height correlation 
and the weak localization effect as a function of mag- 
netic field 6 ^ were confirmed in experiments^*^. Re- 
maining discrepancies between theory and experiments 
regarding the temperature dependence of the width of 
the peak spacing distribution 11 and the peak height dis- 
tributions^ 2 - at low temperatures were explained by the 
inclusion of the exchange interaction term of the univer- 



sal Hamiltonia n 13 i 1 . 

However, not all observed features of the peak spacing 
distribution can be explained by the exchange interac- 
tion alone. At low temperatures, the spacing is given by 
the second-order difference of the ground-state energy 
versus particle number. When only charging energy is 
present, the peak spacing distribution is expected to be 
bimodal because of spin effects. The exchange interaction 
(with realistic values of the exchange coupling constant 
in quantum dots) reduces this bimodality but cannot ex- 
plain its absence in the experiment s 11 ! 15 ' 16 ^ . It is then 
necessary to consider the effect of the fluctuating part of 
the interaction beyond the universal Hamiltonian. 

In the Hartree-Fock-Koopmans approach, the peak 
spacing can be expressed in terms of certain interac- 
tion matrix elements, and sufficiently large fluctuations 
of such matrix elements^ might explain the absence of 
bimodality in the peak spacing distribution. It is there- 
fore of interest to make accurate estimates of interac- 
tion matrix element fluctuations in chaotic dots. These 
fluctuations are determined by single-particle wave func- 
tion correlations. In a diffusive dot, such correlations are 
well understood and lead to an O(Afgx) standard devi- 
ation in the interaction matrix element a 19 ' 20 , where A is 
the mean single-particle level spacing. Peak spacing fluc- 
tuations are also affected by a one-body surface charge 
potential induced by the accumulation of charge on the 
surface of the finite dot^. Matrix element fluctuations 
of the two-body interaction and one-body surface charge 
potential are also important for determining the statis- 
tical scrambling of the Hartree-Fock energy levels and 
wave functions as electrons are added to the do t 21 ' 22 . 

Wave function correlations and interaction matrix ele- 
ments fluctuations in a ballistic dot are less understood. 
In Ref. [23| we used a normalized random wave model 
to obtain analytic expressions for interaction matrix el- 
ement variances and covariances in the regime of large 
Thouless conductance gx for a ballistic two-dimensional 
dot. In such a dot, gx ~ kL, where k is the Fermi wave 



2 



number and L is the linear size of the dot (defined more 
precisely as the square root of the dot's area). Since 
kL ~ y/~N where N is the number of electrons in the dot, 
the kL 3> 1 limit in which the random wave model is ex- 
pected to hold is also the limit of many electrons in the 
dot. In the present work, we systematically investigate 
matrix element fluctuations in real chaotic billiards, for 
30 < kL < 70, corresponding roughly to the parameter 
range relevant for experiments (~ 150 — 800 electrons in 
the dot). We show that fluctuations can be significantly 
enhanced due to dynamical effects, e.g., the variance may 
be enhanced by a factor of 3 or 4. Such enhancement 
can help in explaining the peak spacing distribution mea- 
sured in the chaotic dots of Ref*ii. 

On the other hand, the typical fluctuations of matrix 
elements in chaotic dots cannot explain the even broader 
peak spacing distributions in the experiment of Ref. [l7l 
The small dots used in the latter experiment are probably 
non-chaotic (top gates were used), and this has motivated 
us to study fluctuations beyond the chaotic regime. We 
show that a large (i.e., order of magnitude) enhancement 
of the fluctuations is possible in non-chaotic billiards. 



II. CHAOTIC BILLIARDS 

Here we investigate how dynamical effects modify the 
fluctuations of interaction matrix elements beyond our 
findings in the random wave models. Here and in Sec- 
tions IIIII -M we treat exclusively geometries displaying 
hard chaos. [Systems with stable or marginally stable 
classical trajectories will be considered in Sec. EU] To 
this end, we will use a chaotic system shown in Fig. [T]- a 
modified quarter-stadium billiard geometry^, where the 
quarter-circle has radius R and the straight edge of length 
aR has been replaced by a parabolic bump to eliminate 
bouncing-ball modes. Algebraically, the billiard shape is 
defined by 

< y/R < 1 - s (l - -^jp) > 0<x/R<a 

< y/R < a/I - {x/R-af , a < x/R <a + l, (1) 

where R is the radius of the quarter-circle, and a and s 
are free dimensionless parameters. 



The outline of this paper is as follows. In Sec. [HI we in- 
troduce the modified quarter-stadium billiard as a conve- 
nient model for investigating matrix element fluctuations 
in chaotic systems. In Section lllll we consider matrix el- 
ements of the two-body screened Coulomb interaction, 
and find strong enhancement of the fluctuations in com- 
parison with random wave predictions. Semiclassical cor- 
rections due to bounces from the dot's boundaries lead to 
an increase in the fluctuations, but do not correctly pre- 
dict the scaling with kL in the experimentally relevant 
range. Insight into the underlying mechanism of fluctu- 
ation enhancement is obtained by studying a quantum 
map model, which is described in the Appendix. An im- 
portant conclusion is that the expansion in 1/kL, while 
asymptotically correct, can be problematic in quantify- 
ing matrix element fluctuation in the regime relevant to 
experiments. 

In Sec. IIVI we extend our investigation to one-body 
matrix elements associated with the surface charge po- 
tential, and find similar fluctuation enhancements. Going 
beyond the variance, we examine the full matrix clement 
distributions in Sec. [V] and observe deviations from a 
Gaussian shape that are even stronger than the devia- 
tions found in the random wave models. In Sec. IVII 
we study systems beyond the chaotic regime: billiards 
dominated by marginally-stable bouncing-ball modes and 
billiards with mixed dynamics (i.e., partly regular and 
partly chaotic). Finally, in Sec. IVIll we briefly discuss 
some implications of the present work for the quantita- 
tive understanding of spectral scrambling and peak spac- 
ing statistics for quantum dots in the Coulomb blockade 
regime. 
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FIG. 1: A modified quarter-stadium geometry with parame- 
ters a and s is used to illustrate dynamical effects on matrix 
element fluctuations. In the figure, we set the quarter-circle 
radius R = 1. The random wave contribution to the wave 
function intensity correlator C(ri,r2) is schematically indi- 
cated by a dashed line, and a typical dynamical contribution 
by a dotted line. 

We use a quarter-stadium instead of a full stadium 
shape in order to remove symmetry effects. This system 
has been verified numerically to be fully chaotic for the 
range of parameters used. Variation of the bump size s al- 
lows us to check the sensitivity of the results to details of 
the billiard geometry while maintaining the chaotic char- 
acter of the classical dynamics. Furthermore, by varying 
the parameter a, we can control the degree of classical 
chaos. The degree of chaos can be characterized for ex- 
ample by the Lyapunov exponent A, defined as the rate of 
divergence at long times of generic infinitesimally sepa- 
rated trajectories, \r(t) — r'(t)\ ~ |r — r'|e At as |r — r'| — > 
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and then t -»■ oo. For a = 1.00 and 0.1 < s < 0.2, 
the exponent A takes values 0.69 < XTg < 0.74 (here 
Tb = mL/Tik is a typical time scale associated with one 
bounce in the billiard). When a = 0.25, 0.55 < XT B < 
0.56 in the same range of s, indicating that the system is 
somewhat less chaotic for the smaller value of a. Other 
measures of the degree of chaoticity are possible and may 
be more relevant to the problem of matrix element fluc- 
tuations, as we will argue below. In particular, we may 
consider the rate A* of long-time decay of classical corre- 

lations, f(q,p)f(q(t),p(t)) - f(q,p) ~ e~ A *' as t -> oo, 
where f(q,p) is a typical function defined over the clas- 
sical phase space and the average is over an energy hy- 
persuriace 2 ^. Numerically, we find 0.15 < A*T# < 0.20 
for a = 1 and 0.095 < KT B < 0.13 for a = 0.25, for 
the same range of bump sizes s as above, again indicat- 
ing a less rapid approach to ergodicity in the a = 0.25 
geometry. 

An important consideration in the investigation of dy- 
namical systems, as opposed to random wave models, 
is the presence of boundary conditions. Boundary con- 
ditions lead to Friedel oscillations in the average wave 
function intensity at distances 0(1/ k) from a billiard 
boundary. The effect of such oscillations has recently 
been considered in Refs. H^. The choice of boundary 
conditions, e.g., Neumann or Dirichlct, will also be seen 
to have significant effects on matrix element fluctuations, 
particularly on the fluctuations of one-body matrix ele- 
ments. 

Numerical wave functions for several values of the 
billiard parameters a, s and in various energy ranges 
have been calculated using a variation of the plane wave 
method 27 . At each wave number k, a basis consisting of 
plane waves supplemented by a set of Yq Bessel functions 
centered a fraction of a wavelength outside the boundary 
is used; the size of the basis scales linearly with k. Singu- 
lar value decomposition finds at each k the linear combi- 
nation that minimizes the integrated squared deviation 
along the boundary from the selected boundary condition 
(Dirichlet or Neumann). Finally, minima of this devia- 
tion as a function of k indicate the correct eigenvalues of 
the system. Tests of the method include stability with 
respect to changes in the basis size and comparison of 
the resulting density of states with the Weyl formula. 

Statistics are collected by averaging over an energy 
window. A straightforward estimate shows that such av- 
eraging is sufficient to give good results for matrix ele- 
ment variances, i.e., the ratio of signal to statistical noise 
grows with increasing kL. For all numerical results that 
follow, we use energy windows of constant momentum 
width Ak L = 10, e.g., the data point kL = 30 uses all 
wave functions within the window 25 < kL < 35. The 
Weyl formula for the density of states in two dimensions 
implies that the number of wave functions in such a win- 
dow grows linearly with kL. 



III. TWO-BODY MATRIX ELEMENTS 
A. Fluctuation of diagonal matrix elements v a/ 3 

We first study the variance of the diagonal two-body 
interaction matrix elements v a p = w a /3 ;a( 3, associated 
with a pair of electrons in distinct orbitals a ^ (3 in- 
teracting via the screened Coulomb force. Since the 
screening length of the Coulomb interaction in large 
2D quantum dots is much smaller than the dot size, 
the interaction may be modeled as a contact interaction 
v(r, r') = AV5(r — r'), where V = L 2 is the dot's area, 
and the single-particle mean level spacing A serves to set 
the energy scal o 28 ' 29 . We then have 

v a /3 = AV f dr|V Q (r)| 2 |^(r)| 2 , (2) 
Jv 

where the single-electron wave functions ip obey the usual 
normalization condition J y dr \i[;(r)\ 2 = 1. To leading 
order in 1/gr ~ 1/kL, the variance is then given b y 21 ' 23 : 

where 

C(r, r') = W^WWWW ~ W?W JWW (4) 

is the intensity correlator of a single-electron wave func- 
tion at points r and r'. Assuming C(r, r') is described 
by the normalized random- wave model (i.e., the single- 
electron wave functions are normalized as above with no 
boundary conditions), one obtains 

where (3 = 1, 2 corresponds to the presence or absence of 
time reversal invariance (i.e., the absence or presence of 
an external magnetic field), while b g is a dimensionless 
coefficient that weakly depends on the dot geometry 23 . 

We now evaluate the variance of v a/3 versus kL us- 
ing "exact" (numerically evaluated) real wave functions 
in actual chaotic billiards. Typical results are shown in 
Fig. [21 where we note the large enhancement of the bil- 
liard results over the random wave model (dotted line). 
To understand this enhancement, we compare the ex- 
act numerical results for Sv 2 a with the first term on the 
right hand side of Eq. ([3]), in which C(r, r') is taken to be 
the single-wave-function correlator Cbm(r,r') calculated 
numerically for the appropriate billiard system. The dis- 
crepancy is immediately reduced to a ~ 5 — 10% level, 
which is comparable to the 0((kL)~ 3 ) higher-order cor- 
rection expected and observed in the random wave model. 
Thus, the large enhancement of v a p fluctuations over the 
random wave prediction is not due to higher-order terms 
in Eq. ([3]), but instead can be traced directly to a dynam- 
ical enhancement in the intensity correlator Cbm(r, r') 
over the random- wave correlator. 



4 



5x10"^ 

2x1 0" 2 
1x10" 2 

CO. 

S 

^ 5x1 0- 3 



2x10 



1x10"' 



-3 









2x1 0" 2 






1x10" 2 






30 50 70 



30 



40 



50 



60 



70 



kL 



FIG. 2: The variance of v a/ s versus kL (on a log-linear scale) 
for modified quarter-stadium billiards with Neumann bound- 
ary conditions. The solid line is for a — 0.25, while the dashed 
line is for a = 1.00. In both cases, the results are averaged 
over two values of the bump size: s = 0.1 and 0.2. Dotted line: 
analytic random wave prediction, Eq. ([SJ). Inset: the numer- 
ical result for a = 0.25 with the leading logarithmic term of 
Eq. (0 subtracted (solid line) appears to fall off as (fcL) -1 ' 15 
(dashed line). The analytically expected subleading behavior 
(kL)~ 2 is indicated by a dotted line for comparison. 



We next estimate the dynamical enhancement of the 
intensity correlator (as compared with a random wave 
model) in a semiclassical approach. The random wave 
correlator C rw (r, r') may be interpreted semiclassically 
as arising from straight-line free propagation 20 indicated 
by the dashed line in Fig.[TJ As discussed by Hortikar and 
Srednicki^ - and more recently by Urbina and Richter— , 
additional contributions to the correlator can be asso- 
ciated with trajectories that bounce off the boundary n 
times on their way from r to r', such as the one indicated 
by a dotted line in Fig. [TJ To find these contributions, 
we start from the dynamical correlator for wave function 
amplitudes, which may be written a o 30 ' 31 



■>p*(r)ip(r') 



G*(r,r',E) - G(r', r, E) 
2Trip(E) 



(6) 



where G is the ensemble-averaged part of the retarded 
Green's function G(r,r',E) = ^-Z+u ' and ^ 
is the smooth part of the density of states p(E) = 
Y< a 5 ( E ~ E <*)- Using Eqs. @J|, ©, the dynamical in- 
tensity correlator is given by 



C biU (r, r') = - |G(r', r, E) -G* (r, r' , E) | 2 /4tt 2 f (E) . (7) 

Semiclassically, the smooth density of states is given 
to leading order by the Weyl formula in two dimensions 



p(E) = mL 2 /2irh 2 



while the Green's function is given by the Gutzwiller for- 
mula 3 ^ 



/2 iSj /h—iLij7T /2 



(9) 



The sum in ^ is over classical trajectories j connecting 
r to r' at energy E, Sj is the action along the trajectory 
j, fij is the corresponding Maslov index, and Dj is a clas- 
sical focusing factor that scales as m 2 /pLj (where Lj is 
the length of the trajectory). For the straight-line tra- 
jectory, \Dj\ — m 2 /p\r — r'|. Inserting the semiclassical 
expressions ([8]) and ([9]) into Eq. ([6]), we obtain 



V>*(r)?/>(r') 



1 

V 



J (fc|r-r'|) + Mr,r')(fcLr 1/2 



( 10 ) 

where the Bessel function arises from the straight-line 
path, and h(r,r r ) is a sum over all other trajectories: 



h(r,r>) =J2hj(r,r') 



E 

j 



2pLD } 



Sj (2ns + l)7r" 



(11) 



For typical point pairs (r, r') separated by a distance 
of order L, the function h(r, r') is order unity in kL, 
and the contributions to the correlator from the straight 
line path and from other paths are both 0((fcL) -1 / 2 ). 
For pairs (r,r') separated by a bouncing path of length 
Lj/L < e <C 1, h(r, r') ~ e -1 / 2 . However, the fraction 
of such pairs is 0(e 3 ) and their contribution to the vari- 
ance and other moments of matrix element distributions 
is negligible. 

The intensity correlator in the semiclasssical approxi- 
mation becomes 



C sc (r,r') 



1 2 



J 2 (fc|r-r'|) + / i 2 (r,r')(fc£r 1 



V 2 (3 

+ 2J (fc|r-r'|)Mr,r')(fcL)- 1 / 2 



(12) 



where the first (random wave) term is associated with the 
straight-line path, and the remaining terms constitute 
semiclassical corrections. 

Similarly to the random wave correlato r 23 i 33 ' 34 , 
C sc (r, r') must be corrected to take into account indi- 
vidual wave function normalization. In analogy with 
Refs.llH we have, to leading order in 1/kL, 
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C sc (r,r') = C sc (r,r') + __ 



dr a dr b C sc (r a ,rb) 



v Jv 



-tj [ dr a C sc (r, r a ) - f dr a C sc (r a , r') . 
V Jv v Jv 

Substituting C sc for C in we find 



(13) 



3 /2\ 2 (\nkL + b a ) + b s 



7T U 



(kLf 



O 



A 2 



(kL) 3 



(14) 
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where b sc is a classical constant that in practice must 
be determined numerically by performing the integral in 
Eq. ((3|. As noted above, the random wave and semi- 
classical contributions to C sc (r, r') are of the same order 
except for |r — r'| -C L; it is these short-distance pairs 
that result in a logarithmic enhancement of the random- 
wave term. 

We may easily estimate the dependence of b sc on the 
degree of chaoticity of the dynamical system by invoking 
a diagonal approximation, in which the intensity correla- 
tor C sc (r, r') of Eq. (fT2"j) is averaged over classically small 
regions surrounding r and r'. Noting that Eq. (fTTj) gives 
/z(r,r') as a sum of oscillatory terms with quasi-random 
phases, such averaging leads to 



C s d ( f(r,r') = ^|[^(fc|r-r'|)- 



-E 

j 



/if(r,r')], (15) 



where ■ h 2 (r, r') corresponds to the total classical prob- 
ability of traveling from a neighborhood of r to a neigh- 
borhood of r' via paths j other than the straight-line 
path. Naively, the average semiclassical correction to 
the intensity correlator appears to increase as we in- 
clude longer trajectories. However, let us organize the 
trajectories by number of bounces n or by time t ~ nT B , 
where Tb is a typical time for one bounce in the billiard. 
Trajectories at times t that are significantly longer than 
the classical correlation decay time A" 1 contribute only 
a constant, independent of r and r', to C^ c s (r,r'). This 
is because a classical cloud of trajectories centered near 
r becomes approximately equidistributed over the entire 
billiard when e A ** S> 1, for any initial point r. Such 
position-independent contributions to C~?(r, r') get sub- 
tracted off in the normalization procedure (|13p . Thus, 
the typical size of C sc (r,r') is determined by trajectories 
j having no more than n max w (X*Tg) bounces. 

Furthermore, as a function of t, the number of clas- 
sical trajectories typically grows as e xt , while the focus- 
ing factor for each trajectory j falls off as \Dj\ ~ e~ xt , 
where A is the Lyapunov exponent defined earlier. Thus, 
all n-bounce trajectories combine to form a contribution 
to Eq. (|15[) whose order is roughly n-independent for 
n < "-max- Summing over n up to n max , where n max 
is large, we find 



C s d c g (r 



1 2 



J 2 (fc|r-r'|) + O 



kL 



(16) 

where b\ characterizes the size of the semiclassical con- 
tribution from one-bounce trajectories. Going beyond 
the diagonal approximation is necessary to evaluate prop- 
erly the integral in Eq. © , but the scaling is unaffected. 
Comparing Eqs. ©, (|T1)). and ([To]) , we obtain an es- 
timate for the coefficient b sc in Eq. (fT4)) describing the 
semiclassical correction to the random wave model 



This estimate confirms our intuition that semiclassical 
corrections to the random wave approximation become 
increasingly important as we consider billiards with a 
very long ergodic time A^T 1 . 

Alternatively, the scaling (jTTJ) may be obtained by not- 
ing that when classical correlations persist on a time scale 
A^ 1 that is much longer than the one-bounce time Tb, 
then the effective dimensionless Thouless conductance 
is reduced to gr ~ (\*Ts)kL. Now a typical chaotic 
wave function ip a { r ) may be written as a superposition 
of 0(gx) non-ergodic basis states J?j(r). Since the corre- 
lator r]* (r)Tji(r') for each non-ergodic basis state r]i is of 
order V~ x , we easily see that ^* (r)tjj a (r') takes typical 
values of order V~ 1 g T 1 ^ 2 . The wave function intensity 
correlator C sc (r, r') scales as the square of the amplitude 
correlator, or as V~ 2 g^ 1 for typical pairs (r, r'), yielding 
a lower bound 



{KT B kLf 



(18) 



(17) 



for the integral ©, consistent with Eqs. (H3J) and (JTTJ) . 

For "generic" chaotic systems, the correlation decay 
time A" 1 is of the same order as the one-bounce time 
Tb, and the above scaling arguments for \*T B <C 1 are 
not applicable. Instead, only the first few bounces may 
contribute in practice to the semiclassical correlator, but 
these must be summed up numerically to obtain the semi- 
classical coefficient b sc . This coefficient may in practice 
be quite large even for generic chaotic systems (e.g., the 
modified stadium billiard) and grows as the system be- 
comes less chaotic. 

Qualitatively, the above discussion is consistent with 
our billiard results shown in Fig. [5J as fluctuations are 
observed to be consistently larger for the less chaotic a = 
0.25 billiard, as compared with the a — 1.00 billiard. We 
note that both billiards are "generic" , in the sense that 
they are not fine-tuned to obtain an anomalously long 
time scale A" 1 . We also note that varying the bump size 
s has a very weak effect on the matrix element statistics 
(as long as s is large enough to destroy the bouncing-ball 
modes) and serves instead to provide an estimate of the 
statistical uncertainty in our results. 

For the modified quarter-stadium billiard, we have 
found that adding one-bounce effects to the random wave 
correlator increases the predicted v a /3 variance by ~ 30 - 
40% in the energy range of interest, a significant change 
but not nearly sufficient to explain the full factor of 3 
- 5 enhancement observed in Fig. [3] for the a = 0.25 
billiards (solid lines). Indeed, a close look at the data 
suggests that the numerical results cannot be explained 
fully by semiclassical arguments, no matter how many 
bounces are included in the analysis. The semiclassi- 
cal correction to the variance in Eq. (fl"4"|) is manifestly 
0(l/(kL) 2 ). However, the inset in Fig. [2] clearly shows 
that the dynamical contribution to the variance with kL 
is not consistent with Eq. 114[) but instead appears to fol- 
low a much slower power law ~ l/(fcL)~ 115 . This may 
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FIG. 3: The enhancement of the variance of v a /s (solid line), 
Vac (dashed line) and v a p 1 5 (dotted line) over the correspond- 
ing random wave predictions is shown for a = 0.25 billiards. 
[For v a /3, the random wave prediction is given by Eq. and 
analogous expressions for the other matrix elements may be 
found in Ref. |23| .] In each case, the data is averaged over 
bump sizes s = 0.1 and 0.2. 



FIG. 4: The two-body matrix element variance S for a quan- 
tum map, Eq. (|A6f) in the Appendix, as a function of the 
Hilbert space dimension N. From top to bottom, the three 
solid lines represent data for dominant orbit stability expo- 
nent Ao = 0.25, 0.50, 1.00. The three dashed lines indicate the 
asymptotic 1/N 2 behavior for each case in the semiclassical 
regime of large N. 



be seen also in Fig. [3] (solid line), where the enhance- 
ment over the random wave prediction grows instead of 
diminishing with increasing kL. 

We believe this anomalous behavior results from a 
combination of two related factors: the dynamical en- 
hancement, discussed above, of the b sc coefficient due to 
a finite correlation time scale A^ 1 in an actual dynamical 
system, and the consequent saturation of the l/(kL) 2 
behavior at moderate (< 100) values of kL. As the 
classical system becomes less unstable and the correla- 
tion time A" 1 increases, b sc also increases in accordance 
with Eq. (fTT|) , leading to greatly enhanced matrix ele- 
ment variance at very large values of kL (fT^]) . Because 
the variance is bounded above independent of kL, the 
(fcL)~ 2 growth in the variance necessarily breaks down 
for smaller values of kL. This small-fcL saturation sets 
in at ever larger values of kL as the system becomes less 
unstable and A" 1 becomes larger. 

Alternatively, one may note that the natural expan- 
sion parameter for interaction matrix element fluctua- 
tions in a dynamical system is not (fcL)^ 1 but rather 
the inverse Thouless conductance g^ 1 ~ (X^TskL)^ 1 , 
and the semiclassical contribution with prefactor b sc in 
Eq. (|14p is the leading 0(g^ 2 ) effect in such an expan- 
sion. Terms of third and higher order in g^, 1 , although 
formally subleading and not included in a semiclassical 
calculation, become quantitatively as large as the lead- 
ing 0(g^ 2 ) term when g^ falls below some characteris- 
tic value. Furthermore, if one considers chaotic billiards 
with a long correlation decay time A~ , the importance 
of formally subleading terms in the g^ 1 expansion will 
extend to quite large values of kL. 

The above assertions are explicitly confirmed for a 



quantum map model, described in detail in the Ap- 
pendix, which has scaling behavior analogous to that 
of a two-dimensional billiard, with the number of states 
N = 2n/fi playing the role of semiclassical parameter 
kL = pL/Ti in the billiar d 35 ! 36 . As in the billiard, a free 
parameter in the definition of the map allows for con- 
trol of the classical correlation decay time A" 1 . A key 
difference between the two-dimensional billiard and the 
map model is that the map lacks a logarithmic random 
wave contribution to the variance. We see in Fig. Q] that 
the expected A^ -2 behavior of the variance is observed 
at sufficiently large N, for all three families of quantum 
maps considered. Furthermore, the prefactor multiply- 
ing A^ -2 in each case agrees with that obtained from a 
semiclassical calculation, and as expected this prefactor 
grows with increasing classical correlation time A" 1 (cor- 
responding to a decrease in the chaoticity of the system). 
We also see in Fig. [H that even for a "typical" chaotic 
system (i.e., A*Tb ~ 1), strong deviations from the 1/A 2 
law appear already below N w 80. Such deviations ex- 
tend to even larger N for chaotic systems with slower 
classical correlation decay. This suggests that the large- 
N or large-A:L expansion, though theoretically appealing 
and asymptotically correct, is problematic in describing 
the quantitative behavior of interaction matrix element 
fluctuations for real chaotic systems in the physically in- 
teresting energy range. 

The above numerical calculations were all performed 
in the presence of time reversal symmetry (f3 — 1). From 
Eq. (I14|) we see that when time reversal symmetry is 
broken {[3 — 2), both the random wave contribution to 
the matrix element variance (the term proportional to 
IrtkL + b g ) and the semiclassical contribution (the term 
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proportional to b sc ) are suppressed by the same factor of 
4. Thus, the dynamical enhancement factor for a given 
dot geometry is necessarily /3-independent in the semi- 
classical limit kL ^> 1. However, the saturation effect, 
which tends to suppress the enhancement as kL is re- 
duced, will be less important when = 2, since the vari- 
ance is smaller in this case. Thus, at any finite value of 
kL, the dynamical enhancement in the variance over the 
random wave model will be greater when time reversal 
symmetry is broken, and one may expect enhancements 
somewhat larger than those shown in Fig. [31 This result 
has been confirmed in the quantum map model. 

B. Fluctuation of v aa and v a p-ys 



We have similarly studied the variance Sv^, a of double- 
diagonal interaction matrix elements and the variance 
^ v a0~fS °f off-diagonal interaction matrix elements for ac- 
tual chaotic billiards. Once again, the random wave pre- 
dictions— must be used as the baseline for comparison. 
In Fig. [3j we show the enhancement factor for these ma- 
trix element variances, together with the corresponding 
data for 8v 2 a p discussed previously. 

In the range 30 < kL < 70 most relevant to experi- 
ment, we observe an enhancement in 8v 2 ia over the ran- 
dom wave prediction that is similar to the enhancement 
in Sv^a in the same energy range. In both cases, the en- 
hancement factor continues to grow, instead of approach- 
ing unity, at increasing kL. This latter fact strongly 
suggests that even at kL = 140, we are still far from 
the asymptotic regime of large gx, where matrix element 
fluctuations would be adequately described by a random 
wave picture supplemented by semiclassical corrections. 
The enhancement at large kL is particularly dramatic 
in the case of 5v^ a fluctuations. On the other hand, 
the variance of off-diagonal matrix elements u Q /37<5 is en- 
hanced over the random wave prediction by at most 10%, 
over the entire energy range considered. This is consis- 
tent with the reasonable expectation that dynamical ef- 
fects lead to particularly strong deviations from random 
wave behavior in a modest fraction of the total set of 
single-particle states, such as those associated with par- 
ticularly strong scarring on unstable periodic orbits 37 . 
Such deviations lead to a significant tail in the v aa dis- 
tribution, but have a minimal effect on the distribution 
of off-diagonal matrix elements, since it is unlikely for all 
four wave functions ip a , ipp tp^, and ips to be strongly 
scarred or antiscarred on the same orbit. 

Indeed, inspection of wave functions ip a associated 
with anomalously high double-diagonal matrix elements 
v aa shows that these wave functions have disproportion- 
ately high intensity on average near the dominant hori- 
zontal bounce periodic orbit, which follows the lower edge 
of the billiard in Fig. [TJ We note, however, that asymp- 
totic scar theory in the kL — > oo limit predicts 0(1/ (kL)) 
corrections to the intensity correlation function in posi- 
tion space and only in a region of size 0(1/ '(kL) 1 / 2 ) sur- 



rounding a periodic orbit. Comparing with the integral 
expression @ for the variance, we see that periodic orbits 
asymptotically contribute to the variance only at order 
l/(kL) 3 , compared to the 0(1/ (kL) 2 ) semiclassical effect 
associated with generic (non-periodic) classical trajecto- 
ries (|14[) . Thus, the relative importance of periodic orbit 
effects on matrix element fluctuations is a finite-fc-L (or 
finite-?i) phenomenon, which cannot explain the quan- 
titative scaling behavior of the variance with kL, and 
which is expected to become irrelevant in the asymptotic 
kL —t oo limit. 



C. Matrix element covariance 8v a g8v a ^ 

The normalized random wave model has been shown 
to produce a covariance 8v a p8v ai that is always nega- 
tive, has size ~ A 2 lnkL/(kL) 3 for small uj = Ep — E u , 
and falls off as (oj/E t )~ 2 ~ (SkL)~ 2 for uj ^> E T , where 
Ei is the ballistic Thouless energy^ 3 -. However, in a dif- 
fusive dot, the same matrix element covariance is found 
to be a positive constant oc A 2 / g T (where gx is the dif- 
fusive Thouless conductance) for energy separations uj 
much smaller than the diffusive Thouless energy E c . This 
diffusive covariance falls off for u 3> E c but remains posi- 
tive as long as uj <C Ti/t, where r is the mean free time^. 
An interesting issue is then the sign of the covariance in 
an actual chaotic system. 

First we note the sum rule^ 3 . 

8v a p8v ai = - (Svgp) 2 . (19) 

/3#7 P 

This sum rule is quite general and holds for either a bal- 
listic or a diffusive dot as long as a completeness relation 
is satisfied within an energy window in which the states 
and 7 reside. The average covariance must therefore 
be negative when averaged over all states and 7 within 
such an energy window. The size of the energy window in 
each case must be at least of size h multiplied by the in- 
verse time scale of first recurrences. In a ballistic system 
this implies an energy window of size at least E$ = Ti/Tb, 
where Tb is the one-bounce time. In a diffusive system, 
the completeness relation requires energy scales larger 
than Eq — Ti/t, where r is the mean free time, and thus 
the positive sign of the diffusive covariance at energy sep- 
arations uj <C Ti/t does not contradict the sum rule (|19[) . 

In actual chaotic billiards, it is in principle possible to 
find positive covariance at energy scales ui <§C E , as long 
as the covariance is sufficiently negative for u> ~ Eq to 
produce a negative average covariance over the full en- 
ergy window that is consistent with the sum rule (| 19|) . 
Such positive covariance can result from scars since ijjp 
and ipj will typically be scarred or antiscarred along the 
same orbits when u> — Ep — E 1 is small. The scar contri- 
bution to the covariance for small uj is 0(1/ (kL) 3 ) (i.e., 
of the same order as the scar contribution to the vari- 
ance) and is formally subleading compared with the neg- 
ative 0(hikL/(kL) 3 ) random wave contribution. How- 
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ever, within the range of kL values relevant to experi- 
ments, the scar contribution can dominate and lead to a 
positive covariance for nearby single-particle wave func- 
tions. 



a 
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m/ E n 



FIG. 5: The covariance 8v a fj8v al is computed as a function 
of energy separation uu = Ep — _E 7 for an ensemble of ballis- 
tic discrete-time maps, described in the Appendix, Eqs. (|A9|) 
and (|A10|) . Here Eo = Ti/Tb, where Tb is the one-bounce 
time. The system size N is 128, and A = 0. The dotted line 
indicates the negative average covariance implied by the sum 
rule lH. 



Unfortunately, it is not practical to calculate the ma- 
trix element covariance in a real billiard, since the num- 
ber of wave functions that can be averaged over is not 
sufficient to obtain a signal larger than the statistical 
noise. We instead obtain good statistics for the covari- 
ance in a ballistic discrete map model, introduced pre- 
viously in the discussion of the variance, and described 
in detail in the Appendix. In such discrete maps, the 
matrix element variance or covariance contains no loga- 
rithmic terms. For generic chaotic ballistic systems (i.e., 
Lyapunov time of the same order as the one-step time), 
we find that the covariance is 0(N~ 3 ) ~ 0((kL)~ 3 ) and 
positive for u) -C Eq — Ti/Tb, but becomes negative at 
lu ~ Eo, in contrast with the random wave prediction 
of an always negative covariance. A typical example for 
N = 128 is shown in Fig. \5\ Here discreteness of time 
implies energy periodicity with period 2ttEq — 2t7H/Tb, 
and thus a maximum energy separation lu = nE Q . In 
Fig. the dotted line indicates the negative average co- 
variance over the entire energy window of size 2ttEq, as 
required by the sum rule (fl9)) . 

It is interesting to compare with the covariance in an 
ensemble of two-dimensional diffusive discrete maps^. 
Typical data is shown in Fig. [S] for an ensemble of diffu- 
sive maps on a 32x32 lattice, with Thouless conductance 
c]t = 12 (solid curve) and gr = 24 (dashed curve). The 
theory predicts a variance scaling as 1/g^ and a covari- 
ance scaling as l/<?^, so Sv a f3Sv ai /Sv^ tj3 should scale as 
1/st in the gx —* oo limit. Just as in the ballistic case, 
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FIG. 6: The covariance 8v a p8v al is computed as a function 
of energy separation uo — Ep — E 1 for an ensemble of diffusive 
discrete-time maps on a 32x32 lattice^. The solid curve cor- 
responds to Thouless conductance gr = 12 (E c /Eo — 0.074) 
and the dashed curve corresponds to gr = 24 (E c /Eo = 
0.147). Here Eo — ft/r, where r is the mean free time. The 
value iv = E c , below which the covariance is expected to ap- 
proach a constant positive value, is indicated by a circle in 
each case. The dotted line indicates the negative average co- 
variance implied by the sum rule (|19p . 



the covariance is positive for small separations to and be- 
comes negative when lu ~ E . The average covariance 
over a maximal energy window of size 2ttE is again neg- 
ative, as predicted by the sum rule (fH)|) and indicated by 
a dotted line. 



IV. ONE-BODY MATRIX ELEMENTS 

When an electron is added to the finite dot, charge 
accumulates on the surface and its effect can be de- 
scribed by a one-body potential energy V(r). The di- 
agonal matrix elements of V(r) are given by v a = V aa = 
J v dr \i(j a (r)\ 2 V(r), and the variance of these one-body 
matrix elements may be computed as 



[ [ drdv' V(r)C(r, r')V(r') 



(20) 



Dynamical enhancement of one-body matrix element 
fluctuations may be studied similarly to the analysis 
of two-body matrix element fluctuations presented in 
Sec. IIIII The leading semiclassical contribution to the 
variance is obtained by substituting the normalized semi- 
classical intensity correlator C^ c s [see Eq. ([T5jl ] for C(r, r') 
in Eq. ([2H)l . We immediately obtain 



A 2 



(3 kL 



O 



( A2 



\(kLf 



(21) 



where c g is a geometry-dependent dimensionless coeffi- 
cient arising already in the random wave model 2 ^, while 
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c sc ~ (A*Tb) _1 is associated with the classical dynam- 
ics. We note that the asymptotic power-law behavior of 
the variance is unchanged from the random wave model, 
and the variance is enhanced only by a /cL-independent 
constant. 
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FIG. 7: The variance of the one-body diagonal matrix el- 
ement v a for modified quarter-stadium billiards (a = 0.25; 
averaged over s — 0.1 and s — 0.2) is plotted as a function of 
semiclassical parameter kL. Solid line: Neumann boundary 
conditions. Dashed line: Dirichlet boundary conditions on 
curved boundaries, and Neumann boundary conditions else- 
where. Dotted line: Analytic prediction for the random wave 
model (given by Eq. i|2ip. including only the c 9 term). 



Numerical data for Sv^ in modified quarter-stadium 
billiards is presented in Fig. [71 and compared with ran- 
dom wave results. The ratio of the actual variance to 
the random wave prediction is shown in Fig. [5J Clearly 
this ratio is not constant but rather grows with kL (as 
was also the case with the v a p variance) , indicating once 
again that at kL w 70 we have not yet reached the asymp- 
totic large-fcL regime where semiclassical expressions be- 
come applicable. The same can be observed by compar- 
ing data for Neumann and Dirichlet boundary conditions 
in Fig. [7J Since Dirichlet wave functions decay to zero 
at distances less than 1/fc from a boundary, where the 
surface potential is especially strong, we expect larger 
matrix element fluctuations for the Neumann boundary 
condition data, qualitatively consistent with the results 
in the figure. However, the fraction of points r so close 
to the boundary is 0(1/ kL), while the surface potential 
V(r) is only enhanced by O^kL) 1 / 2 ) there, so the bound- 
ary condition effect is formally subleading. Nevertheless, 
we clearly see from the figure that in the energy range of 
experimental interest, the boundary condition effect is of 
size comparable both to the dynamical enhancement and 
to the baseline random wave prediction for the variance. 
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FIG. 8: Enhancement factor of the v a variance over the ran- 
dom wave prediction is plotted for modified quarter-stadium 
billiards with Neumann boundary conditions, averaged over 
s = 0.1 and 0.2. Solid line: a = 0.25; dashed line: a = 1.00. 



V. MATRIX ELEMENT DISTRIBUTIONS 

Just as was done previously for the random wave 
models, we can go beyond the variance to investigate 
higher moments of the matrix element distribution for 
actual chaotic systems. A typical distribution for diago- 
nal two-body matrix elements v a p in a modified quarter- 
stadium billiard with a = 0.25 and s = 0.1 is shown in 
Fig. O Since the approach to Gaussian behavior is al- 
ready very slow in the case of random waves, it is not 
surprising to find even stronger deviations from a Gaus- 
sian shape for matrix elements in real chaotic systems at 
the same energies. Thus, for modified quarter-stadium 
billiards with a = 1, the skewness 71 of the v a [3 distri- 
bution grows from 1.95 at kL = 70 to 2.72 at kL = 140, 
while the skewness for the same geometry in the random 
wave model drops slightly from 1.21 to 1.09. Similarly, 
the excess kurtosis 72 increases from 8.3 at kL — 70 to 
20.9 at kL = 140, while dropping from 3.7 to 3.3 in the 
random wave model. Similar behavior is obtained for 
other matrix elements. Clearly, the distribution tails are 
very long, and the assumption of Gaussian matrix ele- 
ment distributions is even less justified for real chaotic 
systems than it was in the random wave model. 



VI. BEYOND THE CHAOTIC REGIME 

In this Section we consider fluctuations of matrix el- 
ements in systems that are not fully chaotic. Here no 
universal behavior is expected but we shall see that in 
such systems the variance can be enhanced much more 
than in fully chaotic systems^. We use the modified 
quarter-stadium billiard [see Eq. with s = or 
a < 0. The choice s = corresponds to the original 
Bunimovich stadium, whose quantum fluctuation prop- 
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FIG. 9: The distribution of diagonal interaction matrix ele- 
ments v a /3 is shown for real random waves in a disk22, (dashed 
curve) and for actual eigenstates in a modified quarter- 
stadium billiard geometry with Neumann boundary condi- 
tions (solid curve) at kL = 70. A Gaussian distribution with 
the same mean and variance as the random wave distribution 
is shown as a dotted curve for comparison. 



erties are dominated by the marginally-stable bouncing- 
ball modes, while a < corresponds to a lemon bil- 
liard, which has a classically mixed, or soft chaotic, phase 
space. 



A. Two-body matrix elements 

1. Fluctuation of diagonal matrix elements v a p 

In contrast with the In kL/(kL) 2 falloff in the v a p vari- 
ance predicted for fully chaotic dynamics by Eq. (TH| , in 
the case of regular or mixed dynamics we expect kL- 
independent matrix element fluctuations of order unity. 
To see this explicitly, suppose that the classical phase 
space consists of one regular and one chaotic region, 
with each wave function uniformly distributed over one 
of the two regions. Projecting these regions onto po- 
sition space, let /(r) be the fraction of the energy hy- 
persurface at r that is part of the regular region, i.e., 
the fraction of momentum directions at r that corre- 
spond to stable trajectories. Then the average regu- 
lar wave function has intensity |?/; rog (r)| 2 = V~ 1 f(r)/f 
at position r, while th e average chaotic wave function 
has intensity |^ch(r)| 2 = V^(l - /(r))/(l - /). Here 
/ = y J v dr a /(r a ) is the total fraction of regular points 
in classical phase space, or equivalently the fraction of 
regular quantum eigenstates in the large kL limit. Then, 
starting with the expression ^ for the two-body matrix 
element we find that on average 



whenever a and (3 are both regular states, to be compared 
with the overall average W^p — A for all states a, (3. 
Clearly, v a p is enhanced by a factor of order unity, since 
the two regular states tend to be concentrated in the 
same region of phase space. Similarly, by replacing / with 
1 - /, we obtain enhanced v af3 = A(f 2 -2~J+ 1)/(1 -J) 2 
when both a and (3 are chaotic, and finally, below average 

interaction matrix elements v a [j — A(f 2 — /)/(/ — /) 
are typically obtained when one single-particle state is 
regular and the other chaotic. Combining these results, 
we obtain the lower bound 



Kb > A 2 



P-I 
J-f 



(23) 



,J a(3 



= AV 



dr 

v V 2 



1 f 2 (r) f 2 

J V^ = A^ (22) 



where the quantity in parentheses is a classical system 
property independent of kL. Unless the local regular 
fraction f(r) is a position-independent constant, this 
quantity is nonzero, and the standard deviation is nec- 
essarily of the order of A, i.e. of the same order as the 
average v a p. We note that Eq. (j2"3"j) is a lower bound 
only, as it assumes that each regular or chaotic state is 
uniformly spread over its corresponding phase space re- 
gion. Any intensity fluctuations within the set of regular 
states or within the set of chaotic states will only add to 
the total matrix element variance. 

The fcL-indepcndcncc of the variance can also be in- 
ferred from the following simple argument: regular- like 
quantum behavior is obtained when the ergodic time A" 1 
becomes of the same order as the Heisenberg time irkLTs 
needed to resolve the spectrum. Then the Thouless con- 
ductance gr ~ feLA*Tg is of order unity and Eqs. (fH|) 
and |n]) imply 5v 2 ap ~ A 2 . 

The constant factor in Eq. ([23]) depends not only on the 
regular fraction / in phase space, but equally importantly 

on the relative size ~ / / f 2 of the position-space region 
in which the regular states live (i.e., the participation 
ratio of the regular states). For example, in the extreme 
case where all regular states live in in area V leg and all 
chaotic states live in the complementary area V — Keg, 
we have / = / 2 = V TCg /V, and Sv^ = A 2 , independent 
of the size of Keg- 

Eq. (|2"3"]l predicts very large enhancement, scaling as 
(kL) 2 j InfcL, of the matrix element variance in mixed dy- 
namical systems, over the random wave prediction. Large 
matrix element fluctuations in the presence of soft chaos 
have previously been observed in Ref. |29 . 

The diagonal matrix element variance Sv^p is com- 
puted as a function of kL for two typical mixed phase- 
space quarter-lemon billiards and shown by dashed lines 
in Fig. [101 As expected, no falloff with kL is observed. In 
Fig. Ill) we see that enhancement of an order of magni- 
tude or more over random wave behavior can easily be ob- 
tained for physically interesting values of kL. The most 
dramatic enhancement is observed for the a = —0.25 
quarter-lemon billiard, which is closer to integrability. 

Behavior intermediate between hard chaos and mixed 
chaotic/regular phase space is obtained in the presence 
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FIG. 10: The variance of v a p for a = 0.25, 1.00 quarter- 
stadium billiards (upper and lower solid lines); a = —0.25, 
—0.50 quarter-lemon billiards (upper and lower dashed lines); 
random waves (dotted line). Neumann boundary conditions 
are used for all four billiards. 
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FIG. 11: Enhancement of the v a p variance as compared 
with the random wave prediction for a = 0.25, 1.00 quarter- 
stadium billiards (solid lines); a = —0.25, —0.50 quarter- 
lemon billiards (dashed lines). See Fig. 1101 



of families of marginally stable classical trajectories, such 
as the "bouncing ball" orbits of the stadium billiard. In 
the quarter stadium billiard (s = in Fig. [I]), exceptional 
states associated with such orbits are concentrated in the 
rectangular region of the billiard and constitute a fraction 
~ lf(kL) 1 / 2 of the total set of states^. When a and (3 
are both bouncing ball states, 5v a p = v a f3 — v a p ~ A, just 
as would be the case for regular states concentrated in 
a finite fraction of the available coordinate space. These 
special matrix elements dominate the variance, leading 
to 



K 



kL ' 



(24) 



the random wave prediction. Numerical data for quarter- 
stadium billiards is shown by solid lines in Figs. [TU] and 
ITTI The stronger fluctuations are observed in the less 
chaotic a = 0.25 stadium. 



2. Fluctuation of and Vap-fS 

A calculation analogous to the one resulting in Eq. 
shows that <5v 2 a must also be 0(A 2 ) and fcL-independent 
for a billiard with mixed phase space. In addition, the 
average is enhanced by an 0(1) factor from its ran- 
dom wave value of 3A (/3 = 1) or 2 A ((3 — 2). In the 
stadium billiard, the absence of a stable phase space re- 
gion ensures that bouncing ball states, with 6v aa ~ A 
and frequency ~ lf(kL) 1 / 2 should dominate the double- 
diagonal matrix element variance: 



8v 2 

w v cy.cn 



A 2 



(kL) 1 / 2 



(25) 



and implying an enhancement factor ~ kL/ In kL over 



The billiard results (not shown) are qualitatively consis- 
tent with the above predictions, although statistical noise 
prevents us from extracting a meaningful power law be- 
havior. 

In contrast, fluctuations in the off-diagonal matrix ele- 
ments v a f} 7 s are relatively little affected by bouncing ball 
orbits or even regular phase space regions. This is due 
to the fact that these elements are zero on average, not 
O(A), and thus an increase by an O(l) factor of some 
matrix elements does not necessarily lead to a large vari- 
ance. We may consider an extreme scenario where each 
eigenstate is located in one of two disjoint regions of area 
V/2. Clearly v a /3jS is non- vanishing only when all four 
states are located in the same half of the billiard. In 
such a case, the typical w 2 ^ is enhanced by a factor of 
8 compared with the random wave prediction, ignoring 
logarithms. Because 1/8 of all matrix elements v a ^s are 

nonzero, the variance Sv^ s is nearly unchanged from 
the ergodic case. The above argument generalizes triv- 
ially to an arbitrary number of wave function classes. 
Numerical data in quarter-stadium and quarter-lemon 
billiards (not shown) confirm that Sv^^g is nearly inde- 
pendent of the classical dynamics in the billiard. Higher 
moments of the Sv a /}js distribution are greatly enhanced 
in systems with mixed phase space, and the distribution 
becomes strongly non-Gaussian. 



B. One-body matrix elements 

In a billiard with mixed classical phase space, 
we expect the one-body matrix element v a 
of the surface charge potential V to average 
J v dr V(r)/(r)/ J v dr /(r) = Vf/J for regular states 
where /(r) is the function defined in Section IVI A 11 
and similarly to average (V — V/)/(l — /) for chaotic 
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states. We then obtain a lower bound for the variance 
analogous to Eq. (|2"3")) . 



8vl > 



{vf-vfY 

1-f 



(26) 



which is 0(A 2 ) and independent of kL. Thus, Eq. 
implies an enhancement by a factor ~ kL over the vari- 
ance for fully chaotic billiards given by Eq. (|2"Tj) . The 
absence of a falloff in the variance with increasing kL 
is consistent with our results for quarter-lemon billiards 
(dashed lines) in Fig. [121 




FIG. 12: The variance of v a for a — 0.25, 1.00 quarter- 
stadium billiards (solid lines); a = —0.25, —0.50 lemon bil- 
liards (dashed lines); random waves (dotted line). Neumann 
boundary conditions are used for all four billiards. 

In the quarter-stadium billiard, bouncing-ball states 
with Sv a ~ A will once again dominate the variance 



Svl 



A 2 



(kL) 1 / 2 1 



(27) 



which is a factor ~ (kL) 1 / 2 enhancement over random 
wave behavior. The decay predicted by Eq. (127]) is not 
observed in the numerical data in the experimentally rel- 
evant range 30 < kL < 70 (solid lines in Fig. [12]), sug- 
gesting once again that the energies are not high enough 
for the asymptotic large-fcL scaling laws to be applicable. 
We do find that enhancement by a factor of 5 to 15 of 
the one-body matrix element variance is quite possible in 
the energy range of interest, when the billiard under con- 
sideration exhibits either soft chaos or marginally stable 
orbits in the classical dynamics. 



VII. SUMMARY AND CONCLUSION 

We have studied fluctuations of two-body and one- 
body matrix elements in chaotic billiards as a function of 
a semiclassical parameter kL, and compared them with 



the normalized random wave model predictions. Under- 
standing the quantitative behavior of these fluctuations is 
important for the proper analysis of peak spacing statis- 
tics in the Coulomb blockade regime of weakly coupled 
chaotic quantum dots. 

Dynamical effects, associated with non-random short- 
time behavior in actual chaotic systems, are formally sub- 
leading for two-body matrix elements, and of the same 
order as the random wave prediction for one-body matrix 
elements. In practice, however, we find that these effects 
can easily lead to enhancement by a factor of 3 or 4 of the 
variance in both one-body and two-body matrix elements 
for experimentally relevant values of kL and in reason- 
able hard chaotic geometries. Somewhat larger enhance- 
ment factors are expected when time reversal symmetry 
is broken by a magnetic field. The size of these dynam- 
ical corrections scales in each case as a power of A,: 1 , a 
time scale associated with approach to ergodicity in the 
associated classical dynamics. Random wave behavior is 
recovered in the limit A^T 1 — * 0. In typical geometries, 
dynamical effects on matrix element fluctuations cannot 
be properly computed in a semiclassical approximation, 
as higher-order terms are quantitatively of the same size 
as the semiclassical expression in the kL range of experi- 
mental interest. We have used a quantum map model to 
investigate the approach to semiclassical scaling at very 
large values of kL as well as the saturation of matrix 
element fluctuations at moderate to small values of kL. 

In the case of the interaction matrix element covariance 
for energy levels that are separated by less than the ballis- 
tic Thouless energy, dynamical effects are not only often 
larger than random wave effects, but are also of oppo- 
site sign, leading to an overall covariance that is positive. 
This is in contrast with the random wave model where 
the covariance is always negative. Nevertheless, the sum 
rule (fT^) is preserved due to large negative covariances 
for more widely separated states. We have discussed an 
analogy with similar behavior in diffusive systems. 

Systems with a mixed chaotic-regular phase space or 
with families of marginally stable classical orbits show 
even stronger enhancement of matrix element fluctua- 
tions as compared with the random wave model. We 
discussed the expected asymptotic scaling with kL of the 
matrix clement fluctuations in these cases, and found it 
to be very different from the scaling found in chaotic sys- 
tems. 

Our results strongly indicate that wave function statis- 
tics in actual chaotic single-particle systems, including 
dynamical effects, are needed to make a proper quantita- 
tive comparison between theory (e.g., Hartree-Fock) and 
experiment. A better understanding of single-particle 
wave function correlations is then essential for the calcu- 
lation of observables in an interacting many-electron sys- 
tem such as the peak spacing distribution in the Coulomb 
blockade regime of a quantum dot. Furthermore, these 
correlations need to be understood beyond the naive lead- 
ing order semiclassical approximation, to allow compari- 
son with experiments, which are generally performed at 
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moderate values of the semiclassical parameter kL. 



Acknowledgments 

We acknowledge useful discussions with Y. Gefen, 
Ph. Jacquod, and C. H. Lewenkopf. This work was 
supported in part by the U.S. Department of Energy 
Grants No. DE-FG03-00ER41132 and DE-FG-0291-ER- 
40608 and by the National Science Foundation under 
Grant No. PHY-0545390. We are grateful for the hos- 
pitality of the Institute for Nuclear Theory at the Uni- 
versity of Washington, where this work was completed. 



APPENDIX A: QUANTUM MAP MODEL 



To understand better the anomalously slow decay of 



5v^g and other matrix element fluctuations in realistic 
chaotic systems, we may consider a toy model (perturbed 
cat map^) that displays very similar behavior and for 
which it is easy to collect good statistics at very large 
values of kL. Define a classical map on the torus (q,p) € 

[-7T,7r) X [— 7T, 7r) by 



qt+i = qt + K'(Pt) mod 2ir 
Pt+i = Pt-V(q t +i) mod 2n . 



(Al) 



The above map may be obtained by stroboscopically 
viewing the periodically-kicked Hamiltonian system 



H{q,p,t)=K{p) + J2 S(t-n)V(q), 



(A2) 



n— — oo 



We choose the kick potential to be a perturbation of an 
inverted harmonic oscillator 



V(q) = — — — Acosq — 5(4 cos g — cos2g) 
+ (7(2 sin q — sin 2q) , 



(A3) 



while the kinetic term governing free evolution between 
kicks is 

P 2 

K(p) = — + Acosp + B(4:Cosp - cos2p) . (A4) 

K(p) is even in p to preserve a time-reversal invariance 
(symmetry class (3 — 1). V(q) and K (p) have been chosen 
so that the map has a period-1 orbit at q = p = 0, with 
stability exponent 



An = cosh 1 



I- A. 



(A5) 



where the approximate form holds for An <C 1. Thus, A 
may be varied to change the stability of the shortest orbit, 
whereas the perturbations B and C, which have no effect 
on the linearized behavior around q = p = 0, allow for 



ensemble averaging while keeping the monodromy matrix 
of the central orbit fixed. 

This map may be quantized using standard tech- 
niques--; the position basis is discrete with spacing h 
due to periodicity in momentum. The Hilbert space di- 
mension, N — 2tv/H, plays the role of the semiclassical 
parameter kL = pLjfi in the billiard system. The double 
integral of Eq. @ must be replaced by a double sum 



S = N A 



N 

E [W\^ 



(A6) 



where c is a constant that ensures 



N- 



N 

E [Whfc 



(A7) 



Note that since we are working in one dimension, wc 
must drop the i = j terms to prevent them from dom- 
inating the sum. Our one-dimensional toy model will 
not reproduce the lnkL/(kL) 2 behavior that is associ- 
ated with the short-distance |r — r'| -C L divergence of 
the two-dimensional correlator. Instead, we can think 
of S as the analogue of the two-dimensional integral © 
with the short-distance part subtracted: 



V 2 I I drdr'C 2 (r,r')-- 
v Jv * 



InfcL 



{kL) 2 (kL) 2 



(A8) 

Numerical results for the map are shown in Fig. |4j We 
observe the expected S — b map /N 2 semiclassical behav- 
ior for large N, and the increase of the prefactor 6 m ap 
with decreasing classical stability exponent Ao (see the 
discussion in Section IlII A[) . Furthermore, we note that 
even for the "typical" case Ao = 1, strong deviations 
from the simple power-law behavior appear for N < 50; 
even larger values of ./V are necessary to observe the cor- 
rect power law for smaller Ao- All the curves saturate at 
S w 0.045, leading to the appearance of a slower than 
l/N 2 decay at moderate N values. Thus, it is not sur- 
prising that a weaker than expected dependence on kL 
is observed for moderate kL values in Section IIII Al 

As noted in Ref. [23l the interaction matrix element 
covariance is suppressed relative to the variance by a fac- 
tor ~ kL or N, and the covariance is not a self-averaging 
quantity. To improve the poor ratio of signal to statis- 
tical noise, we may work with a larger ensemble defined 
by 



V(q) = -|- - Acos q + V lnd (q)Q(\q\ - q Q ) (A9) 



and 



K(p) 



V 



+ Acosp + A- rn d(p)e(|(7|-po), (A10) 



where V rn d(q) and K rn d(p) are random functions, -Krnd(p) 
is even to preserve time-reversal symmetry, and is the 
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step function: 0(x) = 1 for x > and otherwise. The we use A — and q = p = tt/2, but very similar 
local dynamics near the periodic orbit at q = p = is behavior is obtained for other values of the parameters, 
unaffected by the ensemble of perturbations. In Fig. [5j 
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